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Abstract. The problem of estimating the eddy difFusivity from Lagrangian observations in the 
presence of measurement error is studied in this paper. We consider a class of incompressible velocity 
fields for which is can be rigorously proved that the small scale dynamics can be parameterised in 
terms of an eddy diffusivity tensor. We show, by means of analysis and numerical experiments, 
that subsampling of the data is necessary for the accurate estimation of the eddy diffusivity. The 
optimal sampling rate depends on the detailed properties of the velocity field. Furthermore, we show 
that averaging over the data only marginally reduces the bias of the estimator due to the multiscale 
structure of the problem, but that it does significantly reduce the effect of observation error. 
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1. Introduction Many phenomena in the physical sciences involve a multitude 
of characteristic temporal and spatial scales. In most cases it is not only impossible 
to study the behavior of the phenomenon at all scales, but it is also unnecessary, since 
usually one is interested in the evolution of a few variables which describe the dynamics 
at large scales. It is, therefore, important to develop systematic methods for deriving 
simplified coarse grained models that capture the essential features of the systems 
at long scales, while accurately parameterising the small scales. In recent years it 
has become clear that the use of data, together with coarse graining procedures, is 
essential for the accurate parameterisation of small scales |GKS041 ICVE06a1 ICVE06bl 
IHKDSOTl IHS08] . The aim of this paper is to study problems of this form in the context 
of transport of passive tracers. 

We are particularly motivated by the challenge of using Lagrangian float data to 
inform the design of subgrid mixing schemes for advected tracers in ocean models. 
The vast amount of Lagrangian float data available (for example, the ARGO project 
has 3000 floats in current operation [ARGO06] ) presents the opportunity to develop 
data-driven model reduction techniques. Lagrangian data are particularly suitable for 
statistical studies of the transport of passively advected substances in the ocean, with 
the simplest statistical description of transport phenomena provided by the average 
concentration of a passive tracer. 

In this paper, we assume that the Lagrangian trajectories are given by the fol- 
lowing stochastic differential equation: 



x = v{x,t) + V2KW. (1.1) 

Here x{t) gM"^ represents the Lagrangian path, v{x,t) is a (prescribed) incompressible 
velocity field, k is the small-scale diffusivity and W{t) denotes standard Brownian 
motion in M''. More sophisticated models have been proposed for oceanographic 
applications, for example |BM02l [BM031 [Pit02l | Wig05[ IGOPR95] . 
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We wish to extract the coarse grained (large length scale and long time scale) 
dynamics of solutions of equation p.ip . For a wide class of velocity fields (de- 
terministic space-time periodic, Gaussian random fields etc.), it is well known 
[MK99, PS08, BLP78] that, at sufficiently long length and time scales, the dynamics 
of (|1.1[) becomes Brownian and can be described by the eddy difFusivity tensor. 
More precisely, it is possible to prove that solutions of (|1.1[1 converge, under the dif- 
fusive rescaling and assuming that the velocity field has zero mean, to an effective 
Brownian motion 



weakly on C([0,T];R''), where W{t) is a standard Brownian motion on R'' and /C 
denotes the eddy (effective) diffusivity tensor. The tensor K, represents the effec- 
tive diffusivity caused by the interaction of molecular diffusion with the transport 
properties of v. 

Consequently, at large length scales and long time scales the dynamics of the 
passive tracer is governed by an equation of the form 



It is quite often the case (in designing subgrid mixing schemes for example) that we 
only wish to calculate the eddy diffusivity, rather than the detailed properties of the 
velocity field v{x,t) at all scales. It is then necessary to estimate the eddy difFusivity of 
a passive tracer from Lagrangian observations. In this paper we address precisely this 
issue: given a Lagrangian trajectory which is consistent with (jl.ip in the presence of 
observation noise, how can we estimate the eddy diffusivity /C? This problem has been 



studied quite extensively over the last few years |BSG02i IF094[ |Fig94| IBSGM098i 
|VGRM04^. 



More generally, we might also want to estimate other coarse grained quantities 
such as the effective drift, or we might want to consider a space dependent eddy diffu- 
sivity. This is a challenging problem in statistical inference: data sampled from (jl.ip 
is only consistent with (|1.3p at sufficiently large scales. In other words, the difficulty 
stems from the fact that the model (|I.3[) used for fitting the data is the wrong model, 
apart from the large scale part of the data. Furthermore, we do not know a priori 
the length and time scales on which the coarse grained model (|1.3p is valid. On the 
other hand, we can perform statistical inference in a fully parametric setting for (jl.3p , 
since only the eddy diffusivity needs to be estimated; statistical inference for (jl.ip 
would require the non-parametric estimation of the velocity field v{x,t) [CDRS09] . 
Parameter estimation for diffusion processes under misspecified or incorrect models 
has been studied in the statistics literature |Kut04l Sec 2.6.1]. 

The problem of parameter estimation for a model that is incompatible with the 
available data at small scales was studied in [PS07[ IPPSOSal IPPSOSbj for a class of 
fast-slow systems of SDKs for which the existence of a coarse grained equation for the 
slow variables can be proved rigorously. In these papers, parameter estimation for the 
averaging problem 




(1.2) 



X = V2K.W. 



(1.3) 




(1.4a) 
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as well as for the homogenization problem 



da; 1 dU , .dV /, r n 

-j^ ^ -Jo{x,y) + Ji[x,y) + ao{x,y) — + ai{x,y) — , (1.5a) 



-jr ^ —9oix,y) + -gi(x,y) + -f3[x,y) — . (1.5b) 
at e e e at 



was studied. In both cases the goal was to fit data obtained from (|1.4al) or (|1.5ap to 
the coarse grained equation 

which describes the dynamics of the slow variable x(t) in the limit as e— >0. In 
the aforementioned papers, it was assumed that the vector field F(X;6) depends on 
a set of parameters 9 that we want to estimate using data taken from either the 
averaging or the homogenization problem. For the homogenization problem it was 
shown in |PPS08a| that the maximum likelihood estimator is asymptotically biased. 
In particular, it is necessary to subsample at an appropriate rate in order to estimate 
the parameters accurately. Similar issues were investigated for the thermal motion of 
a particle in a multiscale potential [PS07| . It was shown that subsampling is necessary 
for the accurate estimation of the drift and diffusion coefficients. 

Related issues have been studied in the field of econometrics. In this context, 
the question is how to accurately estimate the integrated stochastic volatility when 
market microstructure noise {i.e. additive white noise) is present. It was shown 
in |ASMZ05b[ lASMZOSa] that subsampling reduces the bias in the estimator. It was 
also shown that subsampling combined with averaging and an appropriate de-biasing 
step can lead to an accurate and efficient estimator for the integrated stochastic 
volatility. 

In this paper we will study the problem of estimating the eddy diffusivity from 
noisy Lagrangian observations: 



ytj=xt^+0et^, j = l,...7V, 

where {xo,xi, . . . ,xn} is a set of samples from a trajectory consistent with equa- 
tion (jl.ip . etj are independent A/'(0,1) random variables modelling the observation 
error, and measures the strength of the observation error. 

We will consider time-independent spatially-periodic incompressible velocity fields 
as well as spatially-periodic velocity fields that are modulated in time by a time- 
periodic function, or a Gaussian process. In all of these cases, the rescaled trajectory 
converges weakly to a Brownian motion (|1.2p (see |PS08|, Ch. 13]). The eddy diffusiv- 
ity depends in a highly nonlinear way on the properties of the velocity field v{x,t). It 
can be shown (for the class of velocity fields considered in this paper) that the eddy 
diffusivity K. satisfies the upper and lower bounds (we use the notation /C^ — 
where ^ is an arbitrary vector in R'') |AM91| 

K<IC^<-, (1.7) 

K 

for K sufficiently small and some positive constant C. We will consider the physically 
interesting regime 1. 
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As an eddy difFusivity estimator we will use the quadratic variation 

1 

ICn.S = 2]VJ 51 (^"+1 ~ ® i^n+l - Xn) , (1.8) 

where N is the number of observations which we assume to be equidistant, with the 
distance between two subsequent observations being S, and T = NS. It is well known 
[BR80) that, for an SDE of the form we have the convergence result 

lim {x(j+i)T2-^ - XjT2-'^) ® (^(j+i)T2-" - XjT2-N) = 2kIT, a.s. (1.9) 

where / denotes the unit matrix. If we write equation ()1.9|) with Xi replaced by yi, the 
quadratic variation diverges in the limit as A^— s-oo due to the observation error. In 
view of the bounds (|1.7p . it becomes clear that that the estimator ICn^s underestimates 
the value of the eddy diffusivity in this limit. In particular, when the eddy diffusivity 
scales like n^^ the estimator ()1.8p can underestimate the eddy diffusivity by several 
orders of magnitude. 

The above suggests that in order to be able to estimate the eddy diffusivity from 
Lagrangian data, subsampling at an appropriate rate is necessary. However, it is not 
clear a priori what the sampling rate should be. Roughly speaking, we need to look 
at the data at the scale for which the coarse grained description (jl.3|) is valid. The 
estimation of this time scale is a difRcult dynamical question that has been addressed 
only partially |Fan02|, IHP08j . The diffusive time, the time that it takes for the 
Lagrangian particle to reach the asymptotic regime described by a Brownian motion 
with diffusion matrix /C depends crucially on the streamline topology and is related 
to the scaling of the eddy diffusivity with k. Clearly we have two sources of error: 
measurement error, and the error in the estimation of parameters from reduced models 
using data from the full dynamics which we refer to as the multiscale error. The 
multiscale error is precisely due to the fact that the reduced model is incompatible 
with the data at small scales. 

In this paper we study the small k asymptotics of the quadratic variation (jl.Sp . 
We show, by means of rigorous analysis and numerical experiments, that, unless we 
subsample at an appropriate rate, we cannot estimate the eddy diffusivity from the 
quadratic variation, due to the multiscale error. Additionally, we show that for smooth 
time-independent spatially periodic velocity fields, the scaling of the optimal sampling 
rate with k depends on the detailed properties of the velocity field. Our analysis is 
based on standard limit theorems for stochastic processes, together with careful study 
of a Poisson equation posed on the unit torus. 

From the point of view of statistics, it is clearly not optimal to simply ignore 
most of the available data by subsampling. It is natural, therefore, to try to use 
all data through averaging. We experiment with two different types of averaging: 
box averaging (computing the quadratic variation using local averages), and shift 
averaging (which is related to the moving averaging method of statistics). We show 
by means of numerical experiments, that shift averaging significantly reduces the 
effects of observation error, but only marginally reduces the multiscale error. On the 
other hand, box averaging increases the bias of the estimator. 



^We remark, however, that the small scale data that we ignore are highly correlated and it is not 
clear how much additional information they contain about the eddy diffusivity. 
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We emphasize that the setting in which we are working is related to but dif- 
ferent from the problems studied in |PS07[ [PPSOSai TPPSOSb) . In particular, we do 
not assume a priori that we have scale separation and that we know the value of 
the parameter e which measures the degree of scale separation. Rather, the scale 
separation is induced by the dynamics of (|l.ip . The time scale at which the coarse 
grained description is valid is essentially what we need to estimate, since this provides 
us with information about the appropriate sampling rate. For completeness, we also 
consider the rescaled problem (|2.ip below. The rescaled problem and the original one 
are equivalent under space-time rescaling for time-independent velocity fields. How- 
ever, from the point of view of estimating the eddy diffusivity they lead to different 
problems. When using the quadratic variation to estimate the eddy diffusivity in 
equation (|l.ip we actually study the small k limit, whereas for the rescaled problem 
we study the limit of infinite scale separation while keeping k fixed. 

The rest of the paper is organized as follows. In Section[2]we we study the problem 
of estimating the eddy diffusivity using Lagrangian observations from the rescaled 
equation. In Section [3] we study the same problem for the unsealed equation (|l.ip . 
and we also study the effect of observation error on the estimator. In Section [4] 
we develop estimators for the eddy diffusivity which are based on a combination of 
subsampling with averaging and present numerical results for various types of two- 
dimensional velocity fields. Summary and conclusions are presented in Section [5l 
Some technical results are included in the appendices. 

2. The Rescaled Problem 

We consider the equation for the rescaled process 



where we have dropped the superscript e for notational simplicity.. Our goal is to 
estimate the eddy diffusivity using data from (|2.ip . in the parameter regime e<Cl 
and for k fixed. In particular, we want to find how the sampling rate should scale 
with e for the accurate estimation of the eddy diffusivity using the quadratic variation 
estimator. The main result of this section is that, provided that the sampling rate 
is in between the two characteristic time scales 1 and e of the problem, then the 
estimator (|1.8p is asymptotically unbiased, in the limit as e^O. 

We assume that the velocity field is smooth, divergence-free, mean zero and 
1— periodic, i.e. periodic with period 1 in each Cartesian direction. Under these 
assumptions, the solution to (|2.ip converges weakly on C([0,T];R^) to X, as e^O, 
the solution of 



given by 




(2.1) 




dt 



dt 



as e— >0. The eddy diffusivity is given by the formula 




(2.2) 



where the vector field x(z) is the solution of the PDE 



(2.3) 



6 Estimating Eddy DifFusivitios from Noisy Lagrangian Observations 

on T'' with periodic boundary conditions, and where is the generator of the Markov 
process z on T'^: 

- = .(.) + V2^— , (2.4) 

i.e. 

£o = w(2^)-V2 + kA2, 

with periodic boundary conditions. We refer to |PS08) for the derivation of this result. 
Now let /C^ = 5 • K.^ where ^ e K'' arbitrary. From l|2.2p it easily follows that 

K.^^k( |e + V,x«l'd^, (2.5) 

where = X'C- Let /CJ^ g be the quadratic variation along the direction ^: 



N-l 

^b-^E(-»+i--i)'' (2-6) 



n=0 



where a;| = x{n5) ■ ^. 

Our goal is to find how the sampling rate 6 should be chosen so that we can 
estimate the component of the eddy diffusivity (|2.2p along the direction ^ using (|2.6p . 
The following theorem states that the estimator converges in to the eddy diffusivity 
in the limit e— ^0, iV— >oo, with T fixed. 

Theorem 2.1. Let v{z) be a smooth, divergence- free, mean zero, 1-periodic vector 
field and assume that the process z defined in (|2.4p is stationary. Then 

E\IC%,-IC^\^<^ + C{e^S-'' + e^S-^/^ + e^6-^+eS-'/^). (2.7) 
In particular, when 6 = e"', a ^{0,2), we have 

hm limE|/C^_,-/C«|2-0, 

for N6 = T fixed (i.e. iV-e"";. 

Remark 2.1. The scaling of the optimal sampling rate with e, d^e", with a€(0,2) 
appears to be .sharp and it is expected on intuitive grounds, since one would expect 
that the optimal sampling rate should be in between the two characteristic time scales 
of the problem 1 and . 

Remark 2.2. The stationarity assumption on z can be removed .since even when 
z .starts with arbitrary initial conditions its law converges exponentially fast to the 
invariant measure of the process which is the Lebesgue measure on T'^. We refer 
to Wat99^ for the details. 

For the proof of this theorem we will need the following lemma. 
Lemma 2.2. Letv{z) be a smooth, divergence-free, mean zero, 1-periodic vector field 
and assume that the process z defined in (|2.4p is stationary. Then 

In particular, when S — e^^adi (0,2) we have 

lim|E/Ci,_,-/C«|=0. 
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Remark 2.3. Notice that in order for the expectation of the quadratic variation to 
converge to the eddy diffusivity it is not necessary to take the limit N -^oo. Of course, 
in order to keep NS = T fixed we need to take iV~5~^ = K~". 

Proof of Lemma \ 2.2\ With the help of the auxiliary process z — x/e^T'^, equa- 
tion (|2.ip can be rewritten as a system of SDEs: 

dx 1 

— = -v(z) + V2^W. (2.8a) 



-^-y(z) + J-W. (2.8b) 



The generator of the Markov process {x{t), z{t)} is 

£^ = 1 [v{x) • + kA^) + i {v{x) ■ + 2kV^ • V^) + kA^ 
='■ + —Ci+ £2- 
Let x^{^) denote the solution of the Poisson equation 

on T'' with periodic boundary conditions. From standard elliptic PDE theory we have 
that £ C°°{T'^). Hence, we can apply Ito's formula to and use the fact that x^ 
is independent of x to obtain 



dx^ 



« = (^Cox.^ + -Cix^ + C2X^) dt + ^Vyx^ ■ dW 



l-v^[z)dt+^^V ,x^ -dW. 



Consequently: 



1 



f-(n+l)5 An+l)S 

v^{zs)ds = -e{xHzn+i)~xHzn))+V2^ 

^ J n6 J n6 

Thus: 



i+1-4 = -e{xH^n+i)~xHzn))+V2^ / {^zX^+C)-dW 



=:eRn + V2Mn. 
The quadratic variation becomes 



n=0 



Since we have assumed that z(t) is stationary, we have that 

E|M„|2 = .j||e + V,x^||2,(^,)5 = /C«5 (2.9) 
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from which it follows that 



\ n— / 

Furthermore, the maxuTium principle for ehiptic PDEs imphes that 

E|i?„|2<c. 

We use now the above calculations and Cauchy-Schwarz to obtain 

JV-l 

EIC% s ~ -"^^ = ^ E '^^^n + 2V2e]E(i?„M„) 

n=0 

□ 

Proof of Theorem \2.1\ 

From Lemma [2.21 we have that 



E| /Ci,^, - /C« p = E| /C«^ , P - |/C« p + 2/C« (/C« - 

<E|/C«^,|2-|/C«|2 + C(e2ri+e<5-i/2). 

Hence, it is sufficient to estimate the difference i?|/C^ — I^^P- Using the notation 
introduced in the proof of Lemma 12.21 we can write 



N-lN-l 

I^Ll' = ^]^E E (e'^'+2V2ei?„Afn + 2M2) (e^ Rl + 2V2eR,Nh + 2Mf 



n=0 l=Q 
N-lN-1 



i?2 E E^'^«^^'+^' (2.10) 



n=0 1=0 



where 



N-lN-l 



ri=0 1=0 

= :I + II + III + IV + V. 

The uniform bound on a-nd its derivatives, bounds on moments of stochastic inte- 
grals |KS91j and the Cauchy-Schwarz inequality yield the bounds 

EI<Ce'^S-^, EII<Ce^S-^/^, EIII<Ce^S-\ EIV<Ce^6-\ EV<Ced~^/^. 

From the above bounds we deduce that 

ER<C{e*5-^ + e'S-^/^ + e^S-^+e5-^/^). (2.11) 

Now we use bounds on moments of stochastic integrals, together with the fact that 
E{MnMe) = ior nj^l to calculate 

CAT-lAf-l \ N-1 N-1 

^ E E MnM! = ^ E ^(^n) + J^i: Y.nMl)E{Mj) 
n=0 1=0 / n=0 n=0 l^n 
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n=Q i^n 



On the other hand, from Equation (|2.9p we deduce that 

Af-lAf~l 



E 



n=0 e=o 

N-1 



1 "-^ / 1 \ 

— J2J2nM'jEiM!) + o(-). (2.12) 



We combine the above estimates to obtain 

E\IC%,-JC^\'<E\ICif-\Ki\' + C{e'5-^+e6-'/') 



□ 



3. Small K Asymptotics for the Quadratic Variation 

In this section we consider the original problem 

- = .(.) + V2^— . (3.1) 



Our goal is to estimate the eddy diffusivity using data from p.l|) . in the parameter 
regime k<^1. In particular, we want to find how the sampling rate should scale with 
K for the accurate estimation of the eddy diffusivity using the quadratic variation 
estimator. The main result of this section is that in order for the estimator (|1.8p to 
be asymptotically unbiased in the limit as k— >0, it is necessary that the sampling rate 
(as well as the number of observations, and hence the time interval of observation) 
must scale with k in an appropriate way, which depends on the detailed properties of 
the velocity field. In particular, the optimal sampling rate might become unbounded 
in the limit as k — > for flows for which the eddy diffusivity also becomes unbounded 
in this limit. Furthermore, our results are not sharp and detailed analysis is required 
for each particular flow. In contrast with the rescaled problem that was studied in 
the previous section, there doesn't seem to be a simple intuitive argument to explain 
the scaling of the optimal sampling rate with k, since the longest characteristic time 
scale of the problem (the diffusive time scale) needs to be estimated, as a function of 

K. 

As in the previous section we are interested in analyzing the quadratic variation 
along an arbitrary direction ^ and to calculate the optimal sampling rate in order to 
be able to estimate the eddy diffusivity from observations. Let /C^ g be the quadratic 
variation along the direction ^ is given by Equation (|2.6p 

N-l 

'^N.s = ^S^{4+,~4f (3.2) 

n=0 

where x^ = x{n6)-£,. The eddy diffusivity along the direction ^ is given by Equa- 
tion (^3)1 
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where — X' is the unique mean zero solution ol the elliptic PDE 

-{v{z)-V, + kA,)x^^v^ (3.3) 

with periodic boundary conditions on the unit torus. In order to study the small 
K asymptotics of the quadratic variation /C^ g we need information on the small k 
asymptotics of x^, the solution of (|3.3p . From the PDE (|3.3p and Poincare's inequality 
we deduce the bounds 

\\xHl^-<C\\V.xHl^<-. 

K 

The precise asymptotic behavior of x^ in the small k regime depends on the detailed 
properties of the velocity field v{z). This difficult problem has been studied quite 
extensively |CKRZ97[ |BGW89[ IFan02[ IMM93| . In this paper we will assume that the 
solution of the cell problem satisfies the following small-K scaling 

I1x^||lp~||V,x«|Up~k", aehl,0], «;«1, (3.4) 

for p = 2, 4. The notation f ^ n" means that there exists constants C+, C_ so that 

C-K°'<f<C+K°', forK<l. 

Some examples of flows for which the scaling of x^i the solution of (|3.3p . with k is 
known are: 

1. The two-dimensional shear flow v(x) = (0,sin(x)) |MK99i [MM931 . For this 
flow we can solve the Poisson equation explicitly: 

Xi{x,y) = 0, X2{x,y)^~K~'^sm{x) 

and, consequently, for all k>0, 

IIx2||lp~||Vx2||lp~a«-i. 

2. The Taylor-Green flow 

v{x,y) = W^ijjTG {x, y) , (j)TG {x, y) = sin(a;) sin(y). 

In this case it is not possible to solve p.3p . However it is possible to obtain 
sharp estimates on the solution of the Poisson equation: 

||Vx«||l2~^-i/^ k«1 

for all vectors ^gR^. See |Hei03| for details. On the other hand, by the 
maximum principle we have that ||xlU^ ^ ^ i uniformly in k [Fan02] . 

3. The Childress-Soward flow 

v(a;, J/) = V^Vcs (a;,?/), 0CS (a;, y) =sin(a;)sin(?/)-l-Acos(x)cos(?/), 

where Ag[0,1]. This flow interpolates between the Taylor-Green flow (for 
(5 = 0) and a shear flow (for 5 = 1). 
In this case we have that 

\\xHl^-^^-\ WxHl^-I, 
where ^1 = 1/72(1,1), 6 = l/\/2(-l,l). See [SC90l [Fan02] for details. 
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More examples of flows for which the small-K asymptotics of can be calculated will 
be presented in Section ID 

Remark 3.1. Notice that the above scaling leads to 
which is consistent with (|1.7p . since 2a + lG [—1,1]. 

Remark 3.2. Of course, the exponent a in p.4p in general depends in the direction 
^ as well as the U' -space, a = a{^,p). For simplicity we will assume that a is inde- 
pendent of p. The analysis presented below can be easily extended to cover the case 
where a = a{p). 

3.1. Convergence results 

In this section we prove the following. 
Theorem 3.1. Let v{z) be a smooth, divergence-free smooth vector field on T*^. 
Assume that the scaling with p = A holds. Then the following estimate holds 

In particular, if N ^ k'^ with > 4a + 2 and 5 ^ k' with 7 < min(4a + 1, 2a, 8a + 1, ^ + 
i). Then 

limE|/C^ .-/C«p = 0. (3.7) 

K — ^0 ' 

Remark 3.3. Estimate (13. 6p is not sharp. See the examples of the steady and 
modulated in time shear flows in the next section. 

Remark 3.4. Notice that T = NS ^ 00 as k goes to 00, and notice that the sampling 
rate may also have to go to 00 depending on the value of a. This is in constrast to 
the rescaled problem, for which convergence occurs as e — s- with T fixed. 

We first prove the following weak convergence result. 
Lemma 3.2. Letv{z) be a smooth, divergence-free smooth vector field onf^. Assume 
that the scaling \S.4\ with p — 2 holds 

In particular, if6 = K^ wii/i 7<min(2a,4a + l) then 

lini|E/Ci.-/C«|=0. 

Proof. 

We apply Ito's formula to to write the increment of the process as 

l'{n+l)S 

xi+^-xi^V2^ {V,X^+0-dW~ixHzn+i)~xHzn)) 

JnS 

= :V2Mn + Rn, (3.8) 

where 

(M„) = k/ \V,x^+^\'^dz and E(M„)=(5/C«. 

JnS 
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Upon combining (|2.6p and p.Sp and taking the expectation we obtain 

/^W-l 1 N-l 

Ti=0 n=0 

We use now p.4p and (|3.5p to deduce that 

<Ck^"+35-^+Ck^"5-i. (3.9) 

□ 

Proof of Theorem \3.1l 

From Lemma 13.21 we have that 

E/C^ ,5 = /C«+i? (3.10) 

with 

\R\<C{K'^°'+h-^+K^°'6'~^). (3.11) 



We can write 

E|/C^_5-/C«p=E|/C^_5|^-(/C0^-2i?/C«. (3.12) 
We introduce the notation 

I /C^ ,5 1 ^ /2 + //2 + 1 11^ + 2III + 2IIII + 2IIIII 

with 



N6 ^ NS ^ 2N5 ^ " 

n— n— n— 

We use (I2TT2)) to deduce that 
Furthermore, 



Scahng 13.41 together with bounds on moments of stochastic integrals imphcs that 
We conclude that 

E/2<|/C«|2 + c1k4"+2. 
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Consequently 

1 



- ^ N. 
From Assumption p.4p we get 

(E|i?„|P)'/^<CK". 

Now we have 

/ 2 ^-1^-1 \ 

\ n=0 k=0 / 

^ 1^ E E (lE|i?„|^)^/^(E|M„r)V4(E|i?,|4)i/4(E|iV4r)V4 



n=0 fc=0 



Similarly, 



T2p E E ^n^fc 1 
n=0 k=0 / 



We use now the Cauchy-Schwarz inequality to obtain the estimates (we use the fact 
that iV>l) 



E(////)<Ck'*"+^5~\ 
E(/////)<CK4"+5(5-i. 

We use all of the above estimates, together with (I3.12p and estimate p.lip . to obtain 

estimate p.6p . 

□ 

3.2. The Effect of Observation Error 

In this subsection we study the small k asymptotics of the quadratic variation 
in the presence of observation error. More specifically, we assume that the observed 
process (along the direction ^) is 

Y,i^=xl+eel, j = l,...N. (3.13) 

The parameter ^ > measures the strength of the measurement noise which we model 
through a collection of i.i.d Af{Q, 1) random variables . , which are independent from 
the Brownian motion driving the Lagrangian dynamics. Since the two sources of noise 
that appear in the problem are assumed to be independent, the analysis presented in 
this section also applies to Equation (I3.13p . In particular, we have that 
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Fig. 3.1. Figure showing statistics for estimators of the eddy dijfusivity for the shear flow. 
The plots show results for various values of the subsampling interval 5 from (left) the maximum 
likelihood estimator il.SV . (centre) the shift- averaged estimator f^.g[ ), and (right) the box-averaged 
estimator {4.2^ . The plots indicate the mean value of the estimators (circular dots), as well as 
the standard deviation (bars) with statistics computed from 1000 realisations of the Lagrangian 
trajectory. The correct value /C = 5.1, and the value of the small-scale diffusivity ft = 0.1 are both 
indicated as horizontal lines. 



In view of estimate p.9p . we have that 



E 



In particular, \i 5 — tC with 7 < min(2a,4Q + 1,0) then 

lim|E/C^,5-/C«|=0. 



We remark that the exponent 7 is different to the one that appears in the statement 
of Lemma I3.2[ in that it must be negative, irrespective of the scaling of the eddy 
diffusivity with n. 

Similarly, in the presence of measurement error, estimate p.6p has to be modified. 
It becomes 



E 



where R is defined in equation p.lOp and estimated in p. lip . We can then use 
Theorem l3.1l to bound the first term on the right hand side of equation p.l4p . Clearly, 
we require that 5 — > 00 for the additional terms (which are due to the measurement 
error) to vanish. 



C.J. COTTER AND G.A. PAVLIOTIS 



15 



Maximum likelihood Shift averaging Box averaging 

0.30- - 0.30- - 0.30 



0.25 



0.20 



0.15 



0.10 



0.05 



0.00 



A. 



0.25 



0.20 



0.15 



0.10 



0.05 



0.00 



\ 



0.25 



0.20 



0.15 



0.10 



0.05 



0.00 




2 -1 0123 -2 -1 0123 -2 -1 0123 

10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 

a a <5 



Fig. 3.2. Figure showing statistics for estimators of the eddy diffusivity for the periodically- 
modulated shear flow with modulation frequency uj = l. The plots show results for various values of 
the subsampling interval S from (left) the maximum likelihood estimator il.8\ }. (centre) the shift- 
averaged estimator \4-3V , o,i^d (right) the box-averaged estimator \4.SI\j . The plots indicate the mean 
value of the estimators (circular dots), as well as the standard deviation (bars) with statistics com- 
puted from 1000 realisations of the Lagrangian trajectory. The correct value /C = 0.125 (3 d.p.), and 
the value of the small-scale diffusivity k = 0.1 are both indicated as horizontal lines. 



3.3. The Two-Dimensional Shear Flow In this section we present some 
results for a particular class of flows for which we can compute the quadratic variation 
explicitly. The purpose of this is to show that the results obtained in Theorem 13.11 
are not sharp. 

For two-dimensional flows of the form 

v(x,y,t) = {0,f]{t)sm{x)), (3.15) 

where r]{t) can be either a constant, a periodic function or a stochastic process, we 
can calculate explicitly the statistics of the quadratic variation of the Lagrangian 
trajectories jAM90[ [McL98l IMK99] . In the appendix it is shown that for r]{t) = 1, the 
quadratic variation along the direction of the shear is 

where T — N6 and the effective diffusivity is 



16 



Estimating Eddy DifFusivitios from Noisy Lagrangian Observations 



Maximum likelihood Shift averaging Box averaging 

0.20 



0.15 



0.10 



0.05 



0.20 



0.15 



0.10 



0.05 




0.20 



0.15 



0.10 



0.05 



0.00 



2-10123 

10 10 10 10 10 10 



O.OOL 



-2-10123 

10 10 10 10 10 10 
<5 



O.OOL 



-2-10123 

10 10 10 10 10 10 
<5 



Fig. 3.3. Figure showing statistics for estimators of the eddy diffusivity for the O U-modulated 
shear flow with parameters a = l, ct = 0.1. The plots show results for various values of the sub- 
sampling interval 5 from (left) the maximum likelihood estimator il.8\) . (centre) the shift- averaged 
estimator and (right) the box-averaged estimator \4-Sf^ - The plots indicate the mean value of 

the estimators (circular dots), as well as the standard deviation (bars) with statistics computed from 
1000 realisations of the Lagrangian trajectory. The correct value /C = 0.145 (3 d.p.), and the value 
of the small-scale diffusivity k = 0.1 are both indicated as horizontal lines. 



From the above formula we immediately deduce that 

limE[/C7v.5-/C]=0 

provided that 

6 = K-'^-\ (3.17) 
for e > 0, arbitrary. Furthermore, when (|3.17p holds, we have that 

limK-'(E/CA.,,-/C) = -i-^, (3.18) 

the convergence being exponential in k. 

It is also possible to calculate W\JCn .s ^ fO^^ ■ In particular, we have that 

C5 — Vcqh^S"^ +c((5k) I 

K / 
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Fig. 3.4. Figure showing statistics for estimators of the eddy diffusivity for the Taylor- Green 
flow. The plots show results for various values of the subsampling interval 5 from (left) the maximum 
likelihood estimator \1.8\l . (centre) the shift- averaged estimator \4-3^ , and (right) the box-averaged 
estimator {4.2^ . The plots indicate the mean value of the estimators (circular dots), as well as 
the standard deviation (bars) with statistics computed from 1000 realisations of the Lagrangian 
trajectory. The correct value /C = 0.342 (3 d.p.), and the value of the small-scale diffusivity k = 0.1 
are both indicated as horizontal lines. 



di6^ + d:,- + d(iK^5'^ + d{5K)\, (3.19) 

/ 

where the constants {ci, = 1, . . .6} can be calculated explicitly and c{5K),d{5K) 
converge to a constant exponentially quickly in the limit as -\-oo. From the 

above formula we immediately deduce that 

limE|/CAr5-/C|^ = 

K— ►0 

provided that p.lTp holds, together with N ^ k^"^^^ , e>0. Furthermore, under these 
assumptions on 5 and N we have that 

limK""E|/CAr,5-/C|^=const. (3.20) 

This example shows that Theorem 13.11 is not sharp. Some details of the calculation 
of the first two moments of the quadratic variation for the time independent two- 
dimensional shear flow are presented in Appendix [Bl 
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Fig. 3.5. Figure showing statistics for estimators of the eddy diffusivity for the shear flow, 
applied to the rescaled problem with e = 0.1. The plots show results for various values of the sub- 
sampling interval S from (left) the maximum likelihood estimator \1.8\) . (center ) the shift- averaged 
estimator and (right) the box-averaged estimator \4-Sf^ - The plots indicate the mean value 

of the estimators (circular dots), as well as the standard deviation (bars) with statistics computed 
from 1000 realisations of the Lagrangian trajectory. The correct value K = b.l, and the value of the 
small-scale diffusivity k = 0.1 are both indicated as horizontal lines. 



4. Numerical experiments In this section we illustrate the results of the 
previous sections with some numerical experiments, and we investigate some modifi- 
cations to the eddy diffusivity estimator which we shall describe below. The purpose 
of the numerical experiments that we have performed is to investigate the following 
issues: 

1. The performance of the estimator (|1.8p for the eddy diffusivity as a function 
of the sampling rate for flows with different streamline topologies. 

2. Whether an appropriate averaging procedure can reduce the variance of the 
estimator. 

3. The performance of the estimator (|1.8p for the eddy diffusivity as a function 
of the sampling rate for the rescaled problem. 

4. The performance of the estimator (jl.8p in the presence of measurement noise. 
The main conclusions from our numerical experiments can be summarised as follows: 

1. The variance of the estimator as well as the optimal sampling rate depend 
crucially on the streamline topology of the velocity field. 

2. Shift averaging (see below) marginally reduces the variance due to multiscale 
error of the estimator, whereas box averaging (also see below) introduces 
extra bias into the estimator. 

3. There is an optimal sampling rate for the estimator applied to the rescaled 
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Fig. 3.6. Figure showing statistics for estimators of the eddy diffusivity for the periodically- 
modulated shear flow with modulation frequency u) = l, applied to the rescaled problem with e = 0.1. 
The plots show results for various values of the subsampling interval 5 from (left) the maximum 
likelihood estimator il.SV . (centre) the shift- averaged estimator i4-3\ l, and (right) the box-averaged 
estimator The plots indicate the mean value of the estimators (circular dots), as well as 

the standard deviation (bars) with statistics computed from 1000 realisations of the Lagrangian 
trajectory. The correct value /C = 0.125 (3 d.p.), and the value of the small-scale diffusivity k = 0.1 
are both indicated as horizontal lines. 



problem, but even when using the optimal sampling rate the variance of the 
estimator can be very large. 
4. When the data is subject to measurement noise then subsampling is necessary, 
even in the absence of multiscale error. Appropriate averaging can reduce the 
variance due to measurement error. 

4.1. The Estimators 

We are given a time series of Lagrangian observations of length T, sampled at a 
constant rate At. The number of observations is N = T/At. Our goal is to estimate 
the eddy diffusivity using the quadratic variation (|1.8D 

1 

ICn,S = ^ {Xn+l - Xn) (E> {xn+1 - Xn) , (4.1) 

We will consider both the unrescaled (|3.1[) as well as the rescaled problems (|2.1[) . The 
results presented in Sections ^ and ^ suggest that subsampling at an appropriate 
rate is necessary in order to estimate the eddy diffusivity correctly, using Lagrangian 
observations. In the numerical experiments presented in this section we will take the 
sampling rate to scale either with k (for the unrescaled problem) or with e (for the 
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Fig. 3.7. Figure showing statistics for estimators of the eddy diffusivity for the OU-modulated 
shear ftow with parameters q = 1, a = 0.1, applied to the rescaled problem with e = 0.1. The plots 
show results for various values of the subsampling interval 5 from (left) the maximum likelihood 
estimator il.SV . (centre) the shift-averaged estimator f4-3\ l, and (right) the box-averaged estimator 
14- The plots indicate the mean value of the estimators (circular dots), as well as the standard 
deviation (bars) with statistics computed from 1000 realisations of the Lagrangian trajectory. The 
correct value K = 0.145 (3 d.p.), and the value of the small-scale diffusivity k, = 0.1 are both indicated 
as horizontal lines. 



rescaled problem), according to the results presented in Theorems 13.11 and (|2.ip : 

(5~k", or <5~e", 

for some appropriate exponent a. 

Even if we use (|4.ip with 5 chosen optimally, the resulting estimator is clearly not 
optimal since we are using only a very small portion of the available data. Further- 
more, the variance of (jl.Sp with subsampled data can be enormous, in particular when 

1 or e<C 1. One may attempt to reduce the bias and variance in the estimator 
by making use of all the data. In particular, it is reasonable to expect that subsam- 
pling combined with averaging over the data might lead to a more efficient estimator 
of the eddy diffusivity with reduced bias in comparison to the estimator (|4.ip . This 
methodology was applied in [ASMZOSb} I ASMZ05a| in order to estimate the integrated 
stochastic volatility in the presence of market microstructure noise (observation error). 

The most natural way of averaging over the data is by splitting the data into Nb 
bins of size 5 with 5Nb = N and to perform a local averaging over each bin. We use 
the notation 



xi:=x{{n~l)S+{j-l)At), n = 1, . . .TVs, J = 1, • ■ • J, JNb^N, 
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Fig. 3.8. Figure showing statistics for estimators of the eddy diffusivity for the Taylor-Green 
flow, applied to the rescaled problem with e = 0.1. The plots show results for various values of the 
subsampling interval 5 from (left) the maximum likelihood estimator il.8\) . (centre) the shift- averaged 
estimator 14 -Sf , and (right) the box-averaged estimator \4-Sf^ - The plots indicate the mean value of 
the estimators (circular dots), as well as the standard deviation (bars) with statistics computed from 
1000 realisations of the Lagrangian trajectory. The correct value /C = 0.342 (3 d.p.), and the value 
of the small-scale diffusivity k = 0.1 are both indicated as horizontal lines. 



for the j-th observation in the n-th bin. J ^ 6 /At is the number of observations in each 
bin. The maximum hkehhood estimator (jl.Sp is then computed using the averaged 
values 

1 ^ 

■^n = ~^ ^ ^ I 

3 = 1 

leading to the box-averaged estimator: 

= 5^ e' Ui<» - ^hi] « - ^hi] - (4.2) 

«=0 \ j = l 3 = 1 J \ 3 = 1 3 = 1 J 

A second averaging technique, proposed in [ASMZOSal lASMZOSb] to remove the ef- 
fects of market microstructure noise, is to compute a series of estimators, each using 
a different observation from each bin, and then to compute the average. This is the 
shift-averaged estimator: 

^ J ^ Nb-1 

^Ns.s = 7 E ^ E (^"+1 - <) ® - <) ■ (4.3) 

j—1 Tl— 
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Fig. 3.9. Figure showing statistics for estimators of the eddy diffusivity for the shear flow, 
where A/'(0,0.1) observation noise has been added. The plots show results for various values of 
the subsampling interval 5 from (left) the maximum likelihood estimator il.8\ }. (center ) the shift- 
averaged estimator and (right) the box-averaged estimator \4.Sljj . The plots indicate the mean 
value of the estimators (circular dots), as well as the standard deviation (bars) with statistics com- 
puted from 1000 realisations of the Lagrangian trajectory. The correct value /C = 5.1, and the value 
of the small-scale diffusivity k = 0.1 are both indicated as horizontal lines. 



In all of the tests the box-averaged and shift-averaged estimators were obtained using 
values from every single timestep. Throughout this section, we only consider the 
component of the eddy diffusivity along the direction of the shear, since only that 
component is modified by the flow. 

4.2. The Velocity Fields 

The numerical experiments were performed using the following four different ide- 
alized divergence-free velocity fields in two dimensions: 
1. The two-dimensional shear flow: 



v(x) = (0,sin(a;)), 
for which the eddy diffusivity is is |MK99| 



(4.4) 



2. The periodically- modulated two-dimensional shear flow: 



v(x,t) — (0,sin(a;)sin(a.)i)), 



(4.5) 
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Fig. 3.10. Figure showing statistics for estimators of the eddy difjusivity for the periodically- 
modulated shear flow with modulation frequency aj = l, where A/'(0,0.1) observation noise has been 
added. The plots show results for various values of the subsampling interval & from (left) the max- 
imum likelihood estimator (j.gl) . (centre) the shift-averaged estimator {4-3^ , and (right) the box- 
averaged estimator \4.3I^ . The plots indicate the mean value of the estimators (circular dots), as 
well as the standard deviation (bars) with statistics computed from 1000 realisations of the La- 
grangian trajectory. The correct value /C = 0.125 (3 d.p.), and the value of the small-scale difjusivity 
K = 0.1 are both indicated as horizontal lines. 



with uj>0, for which the eddy difFusivity [MBW96] is 

3. The stochastically-modulated two-dimensional shear flow: 

v(x,t) = (0,?7(t)sin(a:)), (4.6) 
where 77 (i) is an Ornstein-Uhlenbeck process obtained from the equation 

and where /3 is a one-dimensional Brownian motion. The eddy diffusivity is 



IC = K- 



2{K-\-a)a 



(4.7) 



The calculation of the eddy difFusivity for this velocity field is presented in 
Appendix [X] 
4. The Taylor-Green flow: 



v{x,t) = V^i^TG (a;, y), V'TG {x,y) = sin(x) sin(y) . 



(4.; 
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Fig. 3.11. Figure showing statistics for estimators of the eddy diffusivity for the OU-modulated 
shear flow with parameters a = l, (7 = 0.1, where A''(0,0.1) observation noise has been added, with 
e = 0.1. The plots show results for various values of the subsampling interval 5 from (left) the 
maximum likelihood estimator (j.gl) . (centre) the shift-averaged estimator \4-3^ , and (right) the box- 
averaged estimator {4.2]j . The plots indicate the mean value of the estimators (circular dots), as well 
as the standard deviation (bars) with statistics computed from 1000 realisations of the Lagrangian 
trajectory. The correct value ^^ = 0.145 (3 d.p.), and the value of the small-scale diffusivity k = 0.1 
are both indicated as horizontal lines. 



There is no closed formula for the eddy diffusivity for this flow, but it is 
well known CS89, SC90l IChi79l IFan02l IKor04j that the eddy diffusivity is 
isotropic and that 



,1/2 



K<1 



with a formula for the prefactor c*. For this case we obtain a numerical 
approximation to the eddy diffusivity K. using the spectral method described 
in |MM93llPi^ . 

We remark that, whereas in the case of the time independent shear flow the eddy dif- 
fusivity becomes singular as 0, in all other examples the eddy diffusivity vanishes 
in the zero molecular diffusion limit. The rate of convergence of /C to is different for 
the velocity flelds ((i?^ and the Taylor-Green flow From Theorem [Owe 

expect that the different scaling of the eddy diffusivity with k should manifest itself 
in the scaling of the optimal subsampling rate with ko 



■^The analysis presented in Section [3] applies only to time-independent velocity fields, but can 
be easily generalized to cover the case of time dependent velocity fields. In fact, for the velocity 
fields l|4.5ll and l|4.6|l we can analyze directly the quadratic variation without appeal to a general 
theory. See Appendix iBl l. 
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Fig. 3.12. Figure showing statistics for estimators of the eddy diffusivity for the Taylor- Green 
flow, where A/'(0,0.1) observation noise has been added. The plots show results for various val- 
ues of the subsampling interval 5 from (left) the maximum likelihood estimator {TTSp, (centre) the 
shift-averaged estimator and (right) the box-averaged estimator The plots indicate the 

mean value of the estimators (circular dots), as well as the standard deviation (bars) with statistics 
computed from 1000 realisations of the Lagrangian trajectory. The correct value Ar = 0.342 (3 d.p.), 
and the value of the small-scale diffusivity k = 0.1 are both indicated as horizontal lines. 



4.3. Results 

Numerical solutions to were obtained for each of these cases using the Euler- 
Maruyama method with a very small timestep to remove the effects of numerical 
discretisation error. The estimator (II. 8p was then computed for each numerical tra- 
jectory and compared with the correct value. In the case of the averaged estimators 
we used all the data in each bin to compute the averages. These calculations were 
repeated for 1000 realisations of the trajectory with different Brownian motions, and 
mean and standard deviations for the estimator values were computed. 

4.3.1. The Unrescaled Process 

Figure ISTTj shows the results of the three estimators applied to the shear flow for 
various values of 6 with an interval width T=1000, from which the number of bins 
Nb = T/5 for the averaged estimators can be computed. As is consistent with equation 
(|1.9p , the maximum likelihood estimator (|1.8p underestimates the eddy diffusivity, and 
converges to the small-scale diffusivity k for small 6. For larger S, the mean value of 
the maximum likelihood estimator approaches the correct value of the eddy diffusivity, 
but the standard deviation of the estimator becomes large, indicating a large variance 
which means that the probability of accurately estimating the correct value is small. 
In comparison, the shift-averaged estimator does not improve the bias by much and 
the variance is only reduced slightly. The box-averaged estimator increases the bias 
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in the estimator in the sense that it substantially underestimates the eddy difFusivity. 

Figure 13.21 shows the same information for the periodically- modulated shear flow 
with modulation frequency uj. The small S limit is again consistent with equation 
(|1.9p , and the mean of the estimator increases to a maximum which is well above the 
correct value, before decreasing again, with increasing standard-deviations for large 
values of S. The shift-averaging again shows very little improvement in either the bias 
or the variance; the box-averaging reduces the mean towards zero in all cases. 

Figure 13.31 shows the same information for the OU-modulated shear flow with 
parameters a= 1, a = 0.1. The results for the maximum likelihood estimator indicate 
an optimum value for S which corresponds with a maximum of the mean, however 
the standard deviation increases monotonically with 5. There is a small improvement 
in the bias and standard deviation for the shift-averaging, and the box-averaging 
produces a mean which is less than the small-scale diffusivity k for all values of 6. 

Figure [3^ shows the same information for the Taylor-Green flow. We observe, as 
is consistent with our theory, that there does seem to be an optimum sampling rate, 
but the variance is large near the optimal rate, similar to the other cases. 

4.3.2. The Rescaled Problem 

We then repeated all of these computations for the rescaled problem (|2.ip with 
e = 0.1. Figures 13.51 13.61 13.71 and 13.81 show the results for the shear flow, the 
periodically-modulated shear flow, the OU-modulated shear flow and the Taylor- 
Green flow respectively. Each of these flows showed that there is an optimal sampling 
rate for which the mean of the maximum likelihood estimator is close to the correct 
value, and that the standard deviation is not too large at this sampling rate, although 
the standard deviation increases for large sampling rates. This illustrates the result of 
theorem 12.11 the mean of the maximum likelihood estimator converges to the correct 
value as e ^ and the variance converges to zero as the subsampling rate 6 converges 
to zero. 

4.3.3. The Effect of Observation Noise In this section we consider the 
combined effect of the multiscale structure and of measurement noise; measurement 
noise is included using equation (|3.13p . The experiments of section 14.3.11 were re- 
peated, with 9 = 0.1. Figures ISTSl 13. 6| 13. 7( and 13.81 show the results for the shear 
flow, the periodically-modulated shear flow, the OU-modulated shear flow and the 
Taylor-Green flow respectively. These results confirm equation (|3.14p in showing that 
the expectation of the estimators tends to infinity as S tends to for non-zero 9. This 
means that it becomes necessary to subsample even if there is no multiscale error. 
The results also show that for 9 = 0.1, the multiscale error dominates the variance of 
the estimator when subsampling is applied. The shift-averaging technique is effec- 
tive at removing the variance due to measurement error, but not the variance due to 
multiscale error. 

5. Conclusions 

The problem of estimating the eddy diffusivity from noisy Lagrangian observa- 
tions was studied in this paper. Apart from the direct relevance of our findings to 
the problem of the accurate parameterisation of the effects of small scales in oceanic 
models, we believe that this work is also a step towards the development of efficient 
methods for data-driven coarse graining. Problems similar to the ones considered in 
this paper have been studied in the context of data assimilation. For example, one 
might fit data from the full dynamics (i.e. the primitive equations) to the quasi- 
geostrophic equation which is a reduced model which is obtained from the full dy- 
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namics after averaging, in the limit as the Rossby number Ro goes to 0. Our resuhs 
suggest that great care has to be taken when fitting data to a reduced model which 
is not compatible with the data at all scales. This is particularly the case when the 
reduced model is obtained through a singular limit such as i?o— >0. 

In this paper, we considered this problem for a class of velocity fields (divergence- 
free, smooth, periodic in space and either steady or modulated in time) for which it 
can be shown rigorously that a parameterisation of the Lagrangian trajectories exists, 
in terms of an eddy diffusivity tensor. For this class of flows, it was shown, by means 
of analysis and numerical experiments, that subsampling is necessary in order to be 
able to estimate the eddy diffusivity from Lagrangian observations. It was also shown 
that the optimal sampling rate depends on the topological properties of the velocity 
field. 

Parameter estimation methods that combine subsampling with averaging of the 
data (defined as shift averaging and box averaging) were also proposed. It was shown 
that shift averaging is very efficient in reducing the effects of observation error, but 
only slightly reduces the variance of the estimator. It appears that the shift-averaging 
technique is only useful for removing measurement error (or microstructure noise in 
the case of econometrics) and not for reducing the multiscale error, as defined in the 
introduction. On the other hand, box averaging leads to a biased estimator, even when 
the optimal sampling rate is used. This should not be surprising, since in the trivial 
case where the velocity field vanishes (i.e. pure Brownian motion with diffusivity 
k), the expectation of the box averaged estimator is k/J where J is the number of 
points per bin. On the other hand, for the same problem, the expectation of the shift 
averaged estimator is n. 

For efficient accurate coarse graining it is necessary to develop estimators which 
can deal with the multiscale error more efficiently. Appropriate averaging over the 
data appears to be an important ingredient of such an estimator. An alternative 
method has been proposed in |CVE06aJ based on the reconstruction of the generator 
of the observed Markov process; methods that combine subsampling and averaging 
with this approach are currently being developed. 

We believe that our conclusions extend to more general types of velocity fields. 
For example, one can carry out the analysis and numerical experiments presented in 
this paper using the class of incompressible Gaussian random velocity fields that were 
considered in [CC99^ . This appears to be a general class of models to consider since 
one can obtain velocity fields with any chosen energy spectrum. The regularity of 
such velocity fields should definitely play an important role in the statistical inference 
procedure. 

Clearly the calculation of the optimal sampling rate from the data is crucial for 
our approach. It appears that frequency domain techniques are more suitable for 
addressing this issue, and this will be investigated in subsequent publications. 
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Appendix A. Derivation of Formula (|4.7p . In this appendix we derive the 
formula for the effective diffusivity for the OU- modulated shear flow (|4.6p . Homoge- 
nization problems for Gaussian incompressible velocity fields that are given in terms 
of an Ornstein-Uhlenbeck process have been considered in [CX971[PSZ07| . The results 
presented in these papers imply that 
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weakly on C([0,T]; 
and 



where W{t) is a standard one-dimensional Brownian motion 



\dx(, 



(1.1) 
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We have used the notation X := (27rT)^ x R, and and p are the unique solutions of 
the equations 

-L(j)^risui{x), / (fipdX^O, (1.2a) 
Jx 



-£*p = 0, / pdX = l. (1.2b) 
Jx 

We have used the notation dX = dxdydrj and £ is the generator of the Markov process 
restricted on X: 

C = r]sm{x)dy + nd^ + ndy — arjdn + ud^ . 

C* denotes the L^(X)-adjoint, i.e. the Fokker-Planck operator. We can easily solve 
equations p.2ap and (|1.2b[) to obtain 



and 



Consequently: 



pdxdydrj^ — e ^ dxdydn, Z = At:'^\ — — 
Z \ a 



(a;,y,??)==— |— ?7sin(a;). 



\U., = j;^.^-'j/^-osix)fpdx 

G 1 



and 



la (k + q;)^ 



%^\\h:p) = T^^Z'' f (sin(x))Vd^ 
1 



2(K + a)2 

Upon inserting the above two formulas in (II. 1|) we obtain (|4.7p . 
Appendix B. The two-dimensional shear flow. 

In this appendix we study in more detail the problem of estimating the eddy 
diffusivity from Lagrangian observations for a class of two-dimensional shear flows. 
Throughout this appendix we only consider the eddy diffusivity along the direction 
of the shear. The flows that we will consider are of the form 

v{x,y,t) = {0,7^{t)f{x)), (2.1) 

where f{x) is a smooth periodic function and r](t) is either a constant, a smooth 
periodic function of time or a stochastic process, e.g. the Onrstein-Uhlenbeck process 
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As it has already been noted in |AM90[ [McL98[ IMKQQ ] , for this class of velocity fields 
the Lagrangian equations can be solved explicitly. In particular, we have that 

y(t)^y{0)+ f Tj{s)f{x{0) + V2^Wi{s))ds + V2^W2{t), (2.2) 
Jo 

where Wi(t) and W2(t) are one dimensional independent Brownian motions. Hence, 
the formula for the quadratic variation becomes 



1 ^-1 / r{n+l)5 \ ^ 

/Cjv,5 = ^X] (y ^ v{s)f{x{0) + V2i^Wiis))ds + V2i^AW2inS)\, (2.3) 

where AW2{nS) = W2{{n+l)S) — W2{n6). Since f{x) is a periodic function, the calcu- 
lation of the statistics of the quadratic variation can be accomplished by calculating 
the statistics of integrals of trigonometric functions of the Brownian motion. This 
calculation can be done by using properties of integrals of symmetric functions, that 
is functions / : [nS, (n-|- 1)^]'^ i-^R for which f{tai,ta2,---taa) = f{ti,t2:---td) for all per- 
mutations cr of (l,2,...rf). In this way, we can calculate the quadratic variation as 
a function of k and 6 in an explicit form. For simplicity we will consider the case 
T]{t) = 1, f{x) = sin(a;) and x(0) = y(0) = 0. The general case can be treated similarly. 
For the velocity field 

v{x,y)^{0,sm{x)) 

We can calculate the expectation of the 22-component of the quadratic variation 
equation (|3.16p . and hence prove (|3.18p . Since Wi{t) and W2{t) are independent, we 
immediately deduce that 

l'(n+l)5 Mn+l)S 

E[(y„+i-2/„)2] = / / E sin(V2^VKi(si))sin(V2^VKi(s2)) d 

In order to calculate the integral on the right hand side of the above equation (which 
we denote by S), we use trigonometric identities together with the formula for the 
expectation of the characteristic function of a Gaussian random variable to obtain 



ael 



(n+l)(5 j-{n+l)5 



nS J nS 

{n+l)S i'{n+l)S 



^iV2^iaiWi{si)+a2Wi{s2)) 



ds2<isi 



\ ^ / / —2K5y^- T a^a^ min(si ,Si ) J j 

= --2_^aia2 ^ ^',^=1 'J ^ " ^'ds2asi, 

aei •^"'^ •^"'5 

where a =(01,02) and / is the index set {— 1, 1} x {— 1, 1} = { — 1, 1}^. The integrand 
is symmetric in si and S2, and, using properties of multiple integrals of symmetric 
functions, we can write the above integral in the form 

-1 /■{n+l)5 /■{n+l)S 

S ^--^aia2 / / e ^^=^ ^^^^<J ^'ds2dsi. 

» JnS Js, 

Evaluating this formula using Maple gives 

S^2kS + - + ^I - ie~^"("+i)* - ie^^™* + ^e-"(4"+i)* + 2e-"* - 2 ) , (2.4) 
K \ 6 2 3 / 
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2k J 




"2k/ 



from which p.l6p follows upon summation. 

We can also calculate E,\ICn,s — , leading to equation ()3.19|) . and hence ()3.20p . 
We have 

/ 1 \ / 1 \2 

(2.5) 

We have already calculated the expectation of the quadratic variation, and it remains 
to compute the second moment. We have 

N N 
n—l m—1 

n—l " n—lm<n ^ 

— —82^ 

We shall separately compute these two types of terms, namely the diagonal terms 5"" 
and the off-diagonal terms 5*2™- 
First we compute S". 



(n+l)(5 > 

sin(\/2Kp^i(s))ds 



\ K 2k^ \ ^ 2 3 / / 

where (|2.4p has been used. For the calculation of S^^ we use trigonometric identi- 
ties, together with the formula for the expectation of the characteristic function of a 
Gaussian random variable to obtain 



1 An+l)S 
" ri,_i Jsk=nS 



aelk=l 



where a= (ai,a2,a3,a4) and / is the indexing set {—1,1}'*. The integrand in this 
multiple integral is a symmetric function, and hence we may write 

o 4 An+1)& An+1)& An+l)S An+1)& 

aG/fe=l J si=nS J S2=si J sa=S2 J S4,=S3 



xG 

This can be computed using Maple: 

" 26880 K4(e«"'5)^®(e'''5)^® 

261 , 1 1 1 1 45 5 

J K~ H 

64 960 (e'«"'5)*V4(e'"')'^ 3360 (e«"'5)*V4(e'"5)^ 16 

49 1 1 1 1 S 



12 K^e'^^ 2400 K4(e'="'5)'^(e'='5)^ 120 ^3 (e'^"'')'' (e'"')'' 
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19 1 1 1 229 1 

24 ^4 ^gKn<5-)4 480 (eKn'5)^^K4gKi5 288 (e'= e'^'' 

, S , 1 1 

-5/4-T-T+3/4- + 



12 /.j3 (■gKn5^4^^5 ^^(e'*"'^)* 1800 ^4^gren5^4(gK(5)4 768 (e''"'')^^K^ 

After summation, all the terms containing exponentials given rise to terms whicli 
converge to a constant divided by S'^N'^ faster than any polynomial power of Sk as 
kS—^oo. 

Next we compute Since n<m, the term inside the sum is 



E[(y"+'-y")'(y'"+'-2/")']=iE 



A7i+l}S /.(n+l)<5 

/ sin(^/2^VKl(s))ds + V2^ / dM^2(s) 

J nS J 7iS 



/ sin(\/2^W^i(s))ds + \/2^ / 

\J mS J n 



= E 



m5 

(n+l)(5 \ ^ / Am+l)S 



(m+l)S ^ 

d 1^2(5) 



sin(\/2KW^i(s))ds / sin(V2KVKi(s))ds 



+2k5 f 2k5 + - + ^ ( - ie-4«(™+i)* _ ig-4«m5 ^ 
\ K 2k^ V " ^ 

2 ^ 

f;g-K(4m+i)i5 ^ 2e^'*'' — 2 

\ K 2k^ V " 2 

|g-^4«+i)5^2e-«^-2j 

+4^252, 

where (|2.4p has been used again. A similar calculation to that for S^i , making use of 
the fact that m>n gives 

<.(n+l)(5 <.(n+l)(5 i / \ , 



r(n+l)S r{n+l)S 1 / 4 \ 

: / / e-2«^'.^=3'^''^^""'^(^-^^)dsids2ds3ds4. 

<y s:^—7n5 J Sd—mS 



The integrand for the two inner integrals, and the integrand for the two outer integrals 
are both symmetric functions, and we obtain 



'-'21 — ~ 



1 /.(ri+l)(5 f.{7i+l)S I ^ \ 
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X / e-2'^^^=3*H"-+25:.<.».".)dsids2ds3ds4. 



This can be computing using Maple: 

anm_ O ^ , 1 ^ 1 1 

'-'21 — ~^ ^ i„ T 
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200 (e"™'^)^ {e^'^f 150 (e'^™'')'' (e''")'' ' 

After the double summation, all the terms containing exponentials given rise 
to terms which converge to a constant multiplied by {N—l)/d'^N'^ faster than any 
polynomial power of Sk as > oo. 

Collecting terms 

E|/C^,.r = ^-^ + ^ + l-^ + '^^ 

NO'^ \ K"^ K'^ K 

+ {di\+d25\+di5^\+di5'^+d^-+dfiK^6^+d{5K) ) , 

iV \ t\j t\j Ki K 
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where the constants {ci,di;i — l,...6} can be read from the above formulas and 
c((5k), c?(5k) converge exponentially fast to a constant in the limit (5k— >+oo. 

Upon computing the remaining terms in equation (|2.5p we notice that all leading 
order terms are cancelled and we end up with 



E\JCn.5-JC\^^ 




which is precisely equation ()3.19|) . 



